/************************************************************************
*************************************************************************
*************************************************************************

simulate_feasible_bivar2.do

Computes whatever stats are needed for a feasible version of the combined
estimator. This basically means bootstrapping the percentiles of the rk stat.

This version varies constraints along 2 dimensions: fraction of firms 
constrained and intensity of constraints. 

This sub-variant imposes the constraint on a fixed set of firms
and is a bit more careful about the intensive margin

		
*************************************************************************
*************************************************************************
************************************************************************/


args jobid sdset

di "ARGS: `jobid'  | `sdset'"

local logon 1

clear
clear matrix
set more off

if mi("`sdset'") {
	local sdset 1
}


//Store current state of the random number generator
global curseed = c(seed)

//Load random seeds
do "do/rngseeds`sdset'.do"

//Load globals
do "constants/simglobals.do"

if $ncores >  1 {
	set seed ${seed`jobid'}
}
	
	
//Set maximum number of iterations
global prevmaxiter = c(maxiter)
set maxiter 1000


if `logon' == 1 {
	cap log close

	cap mkdir "logs"
	local tmp = string(date("`=c(current_date)'", "D M Y"),"%tdCCYYNNDD")
	cap mkdir "logs/`tmp'"
	local tmp2 = clock("`=c(current_time)'","hms")
	log using "logs/`tmp'/`tmp2'simulate_feasible_bivar2`jobid'"
}




//Targets for fraction of firms constrained: Intensive margin
global consttargets 90 45 

//Targets for fraction of firms constrained: Extensive margin
global consttargets_ext90 50 40 30 20 10
global consttargets_ext45 100 80 60 40 20




global dobs 	1



//Name of folder where we'll save the results
global savename sweeps_feasible_bivar2

/************************************************************************
*************************************************************************
Part 1: Helper Programs
************************************************************************
************************************************************************/

//Load data generator
do "do/programs/gendatsyn.do"

//Load GNR program
do "do/programs/gnr_jpe.do"


//Load additional testing programs
do "do/programs/structural_idtests.do"

//Load ar estimation programs
do "do/programs/ar_functions.do"

//Load sharetest program
do "do/programs/sharetest.do"

/************************************************************************
Program: product
Purpose: Computes binary interaction terms
************************************************************************/
cap program drop product
program product
	syntax varlist
	
	gettoken firstvar resvars : varlist	
		
	//Take the product
	foreach var of varlist `varlist' {
		cap gen `firstvar'X`var' = `firstvar'*`var'
	}
	
	//If we're not on the last word, keep going
	if !mi(word("`resvars'",1)) {
		//Continue recursion
		product `resvars'
	}
end







/************************************************************************
Program: truemoments
Purpose: Computes true moments (must tell it the true production function)

Syntax:
name	-	"ces" "quad" (anything else = cobb douglas)
pm		-	Parameter on M (theta)
pk		-	parameter on K
pl		-	parameter on L
mu		-	elasticity of substitution
************************************************************************/
cap program drop truemoments
program truemoments
	syntax name(name=name), pm(real) pk(real) pl(real) [  mu(real 0) sc(real 0) ]
	
	/************************
	Compute marginal products
	***********************/
	
	/****** CES ********/
	if "`name'" == "ces" {
		//For CES they need to pass mu
		
		if mi("`mu'") {
			error "Must pass mu (elasticity of substitution) for a CES production function"
		}
		
		if mi("`sc'") {
			local sc = 1-`pm'
		}
		
		/* There seems to be an error here
		//Marginal product of K and L
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x''*(1-`pm') * Y^(  (1-`pm'-`mu')/(1-`pm')  ) * exp(omega + eps)^( `mu'/(1-`pm') ) / `X'^( 1 - `mu' )
		}
		*/
		
		tempvar I
		gen `I' = `pk'* K^`mu' +   `pl' * L^`mu'
		
		foreach X of varlist K L {
			local x = lower("`X'")
			
			
			
			//This is actually the elasticity times Y/X
			gen MP_`X' =  `sc' *`p`x'' * `X'^`mu' / `I'  * Y/`X'
		}
		
		//Marginal product of M is just the cobb-douglas expression
		gen MP_M = `pm' * Y/M
		
		di "This is CES..."
		
		
	}
	
	/****** Quadratic ********/
	else if "`name'" == "quad" {
		//Marginal product for K and L is just param * Y/X
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x'' * Y/`X'
		}
		
		//Marginal product of M is slightly different
		gen MP_M = exp(omega + eps) * K^`pk' * L^`pl' * (`pm' - 2*M)
	}
	
	/****** Cobb-Douglas ********/
	
	else {
		//Marginal product is just param * Y/X
		foreach X of varlist K L M {
			local x = lower("`X'")
			gen MP_`X' = `p`x'' * Y/`X'
		}	
	}
	
	
	
	
	/************************
	Moments: Elasticities and Sdev of marginal products
	***********************/
	
	//Elasticity = MP_X * X/Y
	foreach X of varlist K L M {
		local x = lower("`X'")
		
		//This may vary by firm
		gen elas`x'_tru 	= MP_`X' * `X'/Y 
		
		//Exclude the top and bottom 1 percent
		bysort t: cumul MP_`X', gen(MPpct)
		
		bysort t: egen sdevMP`x'_tru 	= sd(MP_`X') if MPpct > .01 & MPpct < .99
		
		//Clean up
		drop MP_`X' MPpct
	}

	
end






/************************************************************************
Program: getconstraints
Purpose: Sets a global equal to the proper values of the constraint set

Arguments:
spec	-	The specification to use

************************************************************************/
cap program drop getconstraints
program getconstraints
	args spec
	
	quietly {
	
		
		preserve
		
		local count = 0
		foreach const of numlist -1(.2)20 {
			spec`spec' 2^`const'

			collapse constrained
			
			
			
			local count = `count' + 1
			local frac`count' 	= constrained in 1
			local const`count' 	= `const'
		}

		clear
		set obs `count'
		gen frac 	= .
		gen const 	= .
		forvalues i=1/`count' {
			foreach var of varlist _all {
				replace `var' = ``var'`i'' in `i'
			}
		}
		
		

		global const`spec' ""
		local conlist "${consttargets}"
		foreach target of numlist `conlist' {
			gen tmp = abs(frac-`target'/100)
			sort tmp
			local tmp1 = const in 1
			
			global const`spec' "${const`spec'} `tmp1'"
			drop tmp
		}
		
		restore
		
	}
end




/************************************************************************
Program: gnr_wrap
Purpose: A wrapper for the gnr program (for use in bootstrapping)
************************************************************************/
cap program drop gnr_wrap
program gnr_wrap
	args spec laginst
	
	cap confirm var logtmp_y
	if _rc == 0 {
		//drop y k l m
	}
	else {
		rename (y k l m) logtmp_=
	}
	
	gen share = pricem*M/Y
	//Assume we know when the true process is linear or nonlinear
	if abs(`spec'-6)<.01 {
		//cap gnr y k l m "GMM" "keep" pricem nonlin
		
		if !mi("`laginst'") {
			cap gnr_jpe Y L K M share, `laginst'
		}
		else {
			cap gnr_jpe Y L K M share
		}
	}
	else {
		cap gnr_jpe Y L K M share, linear `laginst'
	} 
	
	if _rc != 0 {
	
		cap gen elasm_gnr 	= .
		cap gen elask_gnr 	= .
		cap gen elasl_gnr 	= .
		cap gen ar_gnr 		= .
	
	}

	drop share

end




/************************************************************************
Program: makebssamp
Purpose: Makes a single bootstrap dataset
************************************************************************/

cap program drop makebssamp
program makebssamp

	/******
	Step 1: Prepare the original dataset for merging
	******/
	
	//Make a household identifier that is numbered consecutively
	egen idcon = group(id)
	
	//Get the max and min value of IDs
	sum idcon
	local idmax = r(max)
	local idmin = r(min)
	
	//Get the max and min values for time
	sum t
	local tmin = r(min)
	local tmax = r(max)
	
	//Number of observations
	egen tag = tag(idcon)
	count if tag
	local nobs = r(N)
	drop tag
	
	
	//Save the source
	tempfile bssource
	save `bssource'
	
	
	/******
	Step 2: Construct framework of BS sample
	******/
	
	//Draw sample
	clear
	set obs `nobs'
	gen idcon = `idmin' + int( runiform()*(`idmax'-`idmin' + 1) )
	
	//The bootstrap household identifier
	gen id_bs = _n
	
	//Add ts
	expand `=`tmax'-`tmin'+1'
	bysort id_bs: gen t = _n + `tmin'-1
	
	tab t
	
	/******
	Step 3: Bring in data from source
	******/
	merge m:1 idcon t using `bssource'
	
	//These are households in the source that were not sampled
	drop if _m == 2
	
	//These are ts in the for which a sampled household is unobserved
	drop if _m == 1
	drop _m
	
	xtset id_bs t

end








/************************************************************************
Program: bspctiles
Purpose: Bootstraps the percentiles of the rk-statistic (and the vcov of the Lambda stat)
************************************************************************/

cap program drop bspctiles
program bspctiles
	//Number of reps, stub of the estimator (gnr or ar_dem), specification,
	//and a specifier for whether we should save or not
	args nreps spec savepath
	
	/******
	Step 0: Estimate true parameters
	******/
	preserve
	sort id t
	ar_dem
	local rhoest = _b[/rho]
	foreach inst of varlist k-mXm {
		gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
		gen INST_`inst' = l.`inst'
	}
	
	gen yrho2 = y - `rhoest'*l.y
	
	foreach endvar of varlist ENDO2_* {
		reg `endvar' INST_* l2.m l2.k l2.kXm
		estimates store `endvar'
	}
	
	ranktest_hacked (ENDO2_*) (INST_* l2.m l2.k l2.kXm), cluster(id) wald fullrank
	local nel = rowsof(r(lab))
	
	local listel
	forvalues i=1/`nel' {
		local listel `listel' lab`i'
	}
	
	
	//ivreg2 yrho2 (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id)
	restore
	//gen widstat_iv = e(widstat)

	/******
	Step 1: Prepare a dataset for the bootstrapped estimates
	******/
	cap postclose bsests
	
	tempfile bsests
	postfile bsests widnull `listel' using `bsests'
	
	
	/******
	Step 2: Prepare the original dataset for bootstrapping
	******/
	tempfile original
	save `original'
	
	
	
	
	/******
	Step 3: Bootstrap
	******/
	forvalues b=1/`nreps' {
		preserve
		makebssamp
		drop id idcon
		sort id_bs t
		
		//Preliminary construction
		ar_dem
		local rhoest = _b[/rho]
		foreach inst of varlist k-mXm {
			gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
			gen INST_`inst' = l.`inst'
		}
		
		//Rank test, etc.
		ranktest_hacked (ENDO2_*) (INST_* l2.m l2.k l2.kXm), cluster(id_bs) wald fullrank
		matrix l = r(lab)
		local passlist
		forvalues i=1/`nel' {
			local passlist `passlist' ( `=el("l",`i',1)' )
		}
		
		//Nullify
		if "${test_fullrank}" != "fullrank" {
			//Null: zero rank
			foreach endvar of varlist ENDO2_* {
				estimates restore `endvar'
				predict tmp
				replace `endvar' = `endvar' - tmp
				drop tmp
			}
		}
		else {
			//Null: K-1 rank
			foreach endvar of varlist ENDO2_m {
				estimates restore `endvar'
				predict tmp
				replace `endvar' = `endvar' - tmp
				drop tmp
			}		
		}
		
		gen yrho2 = y - `rhoest'*l.y
		ivreg2 yrho2 (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id_bs)
		local widnull = e(widstat)
		

		post bsests (`widnull') `passlist'
		restore
	}
	
	postclose bsests
	
	/******
	Step 4: Compute standard errors
	******/
	
	use `bsests', clear
	
	if !mi("`savepath'") {
		save "`savepath'", replace
	}
	
	local count 0
	foreach p of numlist 10 5 1 .5 {
		local ++count
		egen tmp = pctile(widnull), p(`=100-`p'')
		sum tmp
		local p`count' = r(mean)
		drop tmp
	}
	
	cap matrix drop labv
	corr lab*, cov
	matrix labv = r(C)
	

	
	use `original', clear
	
	foreach p of numlist 1/`count' {
		gen widstatp_`p' = `p`p''
	}
	
	matrix veclabv = vec(labv)'
	svmat veclabv


end



/************************************************************************
*************************************************************************
Part 2: Scenarios
************************************************************************
************************************************************************/


/************************************************************************
Fixed parameters
************************************************************************/

//List of constraints (in log base 2)
//local conlist -3(2)7

//Parameters of the production function
global pk 		= .11
global pl 		= .28
global pm 		= .67


//global quad 	= 10


//CES, almost cobb-douglas
global mu 		= .2   //Implies an elasticity of substitution equal to 1.25
global cespk 	= .064
global cespl 	= .33
global sc 		= .39

//CES, not very cobb-douglas
global mu2		= .8   //Implies an elasticity of substitution equal to 5
global cespk2	= .015
global cespl2	= .95
global sc2		= .39

////Standard deviations

//Anticipated shock
global sigeta1 	= .0925324 
global sigeta2 	= ${sigeta1} * 2

//Unanticipated shock
//global sigeps 	= 1.07
global sigeps 	= .2542193



//Number of firms
global nfirmnorm 	2613
global nfirmsmall 	1000


//Number of years
global nyearnorm 	7
global nyearshort 	4

//Persistance of productivity
global rhocal .7709392



/*********************
Production Functions
**********************/

/*
The DGP program can handle production functions of the form

X(K,L)*M^pm

OR

X(K,L)*(pm1*M - M^2)

where X(K,L) is arbitrary. This works because in both cases optimal M depends on
only X(K,L), but not K and L separately.
*/

global XKL1 "gen X_kl = K^${pk} * L^${pl}"
global XKL2 "gen X_kl = (  ${cespk}*K^${mu} + ${cespl}*L^${mu} ) ^ ( ${sc}/${mu} )"
global XKL3 "gen X_kl = (  ${cespk2}*K^${mu2} + ${cespl2}*L^${mu2} ) ^ ( ${sc2}/${mu2} )"


/*********************
Productivity processes
**********************/
global nlomproc "86.30933 - 7.081727 - 33.35803*(l.omega + 7.081727) + 4.568181*(l.omega + 7.081727)^2 - .2030311*(l.omega + 7.081727)^3 "



/************************************************************************
Normal Specifications
************************************************************************/
cap program drop spec1

//Baseline
program spec1
	args const fracobs
	

	
	gendatsyn  =`const' , klfunc("${XKL3}") rho(${rhocal}) n(${nfirmnorm}) t(${nyearnorm}) pm(${pm}) sigeps(${sigeps}) sigeta(${sigeta1}) 
	
	truemoments ces, pm(${pm}) pk(${cespk2}) pl(${cespl2}) mu(${mu2}) sc(${sc2})
end


/************************************************************************
Bivariate Specifications
************************************************************************/
//Drop programs in memory

//Arguments to the DGP: klfunc mfunc constraintmult rho N T pm sigeps sigeta pricemsd


//Baseline
cap program drop bivarspec1
program bivarspec1
	args const fraccon
	
	local ncon = round(`fraccon'/100*${nfirmnorm})
	local nunc = ${nfirmnorm} - `ncon'
	
	if `nunc' > 0 {
		gendatsyn  =`const' , klfunc("${XKL3}") rho(${rhocal}) n(`ncon') t(${nyearnorm}) pm(${pm}) sigeps(${sigeps}) sigeta(${sigeta1}) 
		tempfile con
		save `con'
		
		gendatsyn  = 2^20 , klfunc("${XKL3}") rho(${rhocal}) n(`nunc') t(${nyearnorm}) pm(${pm}) sigeps(${sigeps}) sigeta(${sigeta1}) 
		replace id = id+`ncon'
		append using `con'
	}
	else {
		gendatsyn  =`const' , klfunc("${XKL3}") rho(${rhocal}) n(`ncon') t(${nyearnorm}) pm(${pm}) sigeps(${sigeps}) sigeta(${sigeta1}) 		
	}
	
	
	truemoments ces, pm(${pm}) pk(${cespk2}) pl(${cespl2}) mu(${mu2}) sc(${sc2})
end



/************************************************************************
*************************************************************************
Part 3: Sweep
************************************************************************
************************************************************************/
local nreps = ${B} 

cap mkdir "dta/simulations"
cap mkdir "dta/simulations/${savename}"

if ${dobs} == 1 {
	cap mkdir "dta/simulations/${savename}_bs"
}

if ${ncores} > 1 {
	cap mkdir "dta/simulations/${savename}/`jobid'_`sdset'"
	cap mkdir "dta/simulations/${savename}_bs/`jobid'_`sdset'"
	
	local savefol "dta/simulations/${savename}/`jobid'_`sdset'"
	local bsfol "dta/simulations/${savename}_bs/`jobid'_`sdset'"
	
	di "*******************"
	di "I am job `jobid'"
	di "*******************" _n _n _n
}
else {
	local savefol "dta/simulations/${savename}"
	local bsfol "dta/simulations/${savename}_bs"
}



/*
Estimate time to complete
How many operations?
*/
//local speclist 1 4 5 6
local speclist 1 
local nspecs = wordcount("`speclist'")

local totops = `nspecs' * wordcount("${consttargets}") * wordcount("${consttargets_ext90}") * `nreps'
local opscom = 0
local esttime "???"



//Loop through specifications
foreach spec of numlist `speclist' {

	cap mkdir "`savefol'/spec`spec'"
	if ${dobs} == 1 {
		cap mkdir "`bsfol'/spec`spec'"
		local bshere "`bsfol'/spec`spec'"
	}
	
	di _n "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"

	
	//Get the constraints for this specification
	getconstraints `spec'
	local conlist "${const`spec'}"
	
	local ind = 0
	foreach fc of numlist ${consttargets} {
		quietly: local ++ind
		local const: word `ind' of `conlist'
		
		foreach extcon in ${consttargets_ext`fc'} {
			
		
			di "-----------------------------------------------------------------------------------------"
			
			
			
			forvalues b=1/`nreps' {
				//Version of the constraint to display
				di "On Specification `spec'" _col(25) "Constraint = (Ext=`extcon',Int=`fc')" _col(70) "Rep = `b'" _col(85) "Est. Time: `esttime' min"
			
				quietly {
					timer on 1
			
					/********
					Generate data and compute true moments
					********/
					bivarspec`spec' 2^`const' `extcon'
					order elas* sdevMP*
					
					
					
			
					
					/*************************************
					Step 1: Apply Methods
					*************************************/
					
					//From here it's almost identical to our earlier approaches,
					//except we're only using the specifications that seem valuable
					
					//Also, we'll do both with and without the lag assumption
					

					
					//Construct second-order terms
					order k l m, last
					product k l m
				
					/*********
					Share test
					*********/
					gen logshare = log( pricem*M/Y )
					sort id t
					//sharetest logshare k l m, test( l2.( c.k##c.m ) ) lagtest deg(4)
					sharetest logshare k l m, test( l2.( c.k##c.m ) ) lagtest deg(2) lagdeg(1)
					
					local stub idtest_share
					gen `stub'_pr2	= r(r2)
					gen `stub'_F	= r(F)
					gen `stub'_pval	= r(pval)
					
					
					
					sharetest logshare k l m, test( l2.( c.k##c.m ) ) mmerr deg(2)
					
					local stub idtest_mmshare
					gen `stub'_pr2	= r(r2)
					gen `stub'_F	= r(F)
					gen `stub'_pval	= r(pval)
					
					
					/********
					WIDstat test
					********/
					widstat_iv_old
					widstat_iv

			
					
					/********
					Get estimates
					********/
					
					foreach estimator in ar_dem ${arselect} gnr {
						sort id t
						`estimator'_wrap `spec' 
						
						/*
						if ${dobs} == 1 {
							local fc : word `ind' of $consttargets
							bserr 50 `estimator' `spec' "`bsfol'/spec`spec'/const`fc'_`estimator'_rep`b'.dta"
						}*/
					}
					
						
					
					
					//Compute standard deviation of marginal products
					//foreach method in ar_std ar_dem gnr {
					foreach method in ar_dem gnr {
						foreach var of varlist K L M {
							local sm = lower("`var'")
							gen MP_`sm' = elas`sm'_`method' * Y/`var'
							
							//Exclude the top and bottom 1 percent
							bysort t: cumul MP_`sm', gen(MPpct)
							
							bysort t: egen sdevMP`sm'_`method' = sd(MP_`sm') if MPpct > .01 & MPpct < .99
							drop MP_`sm' MPpct
						}
					}
					
					
				
					/*************************************
					Step 2: Run Tests of Scalar Unobservable
					*************************************/
					//In reality you'd be running these tests first, but they
					//tend to conjure up a bunch of new variables that may create
					//difficulties
					
					sort id t
					cap drop k l
					rename logtmp_* *
					order k l m kXk kXl kXm lXl lXm mXm, last
					
					
					/********
					Form observation
					********/
					//keep elas* sdevMP* test* widstat widpval constrained 
					if ${dobs} == 1 {
						//local keepse se*
						bspctiles 100 1 "`bshere'/`extcon'_`fc'_`b'.dta"
					}
		
					
					//keep elas* ar_* `keepse' `keepsel' sdevMP* test* constrained jacmat* idtest*  hybrid_test*
					keep elas* ar_*  sdevMP*  constrained idtest*  dofadj widstat_iv* widstatp_* labest* veclabv*
					
					
					collapse _all
					
					tempfile `b'
					save ``b''
					
					
					/********
					Estimate time remaining
					********/
					timer off 1
					timer list
					local opscom 	= `opscom'+1
					local timesofar = r(t1)
					local avgtime	= `timesofar'/`opscom'
					local timerem 	= `avgtime' * (`totops'-`opscom') / 60
					local esttime 	= string(`timerem', "%9.2f" )
					
				}
				
			
			}
			
			quietly {
				use `1'
				forvalues k=2/`nreps' {
					append using ``k''
				}
				
				
				
				//if there's a negative, change it to an underscore
				gen const 	= `const'
				gen fraccon	= `fc'
				

				save "`savefol'/spec`spec'/const`extcon'_`fc'.dta", replace
				
				
			}
		}
	}
	
}


//Restore default max number of iterations for GMM
set maxiter ${prevmaxiter}


if `logon' == 1 {
	cap log close
}


